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ABSTRACT 


This  paper  ia  a  survey  of  finite  difference  techniques  for  solving  the 
incompressible  Navier-Stokes  equations.  We  consider  what  features  are 
important  for  schemes  in  order  to  be  useful  in  scientific  and  engineering 
applications.  In  so  far  as  possible  we  try  to  point  out  the  strengths  and 
weaknesses  of  various  methods. 
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SIGNIFICANCE  AND  EXPLANATION 


\ 


V 

The  incompressible  Navier-Stokes  equations  describe  the  flow  of  fluids 
such  as  water,  air  at  low  speeds,  and  many  other  common  liquids  and  gases. 
Efficient  and  accurate  numerical  methods  for  their  solution  are  therefore  of 
great  importance  for  a  wide  variety  of  engineering  and  scientific 
applications.  In  this  paper,  several  existing  finite  difference  methods  are 
examined  and  suggestions  for  further  research  are  given. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MFC,  and  not  with  the  author  of  this  report. 


FINITE  DIFFERENCE  METHODS  FOR  TIB  INCOMPRESSIBLE 
NAVIER-STOKES  EQUATIONS  -  A  SURVEY 

John  C.  Strikwerda 

1.  Introduction 

The  incompressible  Navier-8tokes  aquations  doacriba  the  motion  of  fluids  such  as 
“•ter,  air  at  low  speeds,  and  sany  other  common  liquids  and  gases.  As  such  they  are  of 
widespread  scientific  and  engineering  aignificance  and  numerical  solution  procedures 
equations  are  of  great  importance.  In  this  paper  we  will  examine  common  numerical  methoda 
for  solving  the  incompressible  Navier-stokes  equations  and  discuss  various  strengths  and 
weaknesses  of  the  methods.  We  will  particularly  concern  ourselves  with  finite  difference 
methods  although  much  of  this  paper  will  apply  to  general  numerical  methods. 

The  incompressible  Navier-stokes  equations  are 

(1.1)  ufc+  (u  •  $)u)  *  $  p  -  V2u  -  f(x,t) 

$  •  u  -  0 

where  u(x,t)  ia  the  velocity  p(x,t)  is  the  pressure,  and  $(x,t)  repreaents  forces  on 
the  fluid.  The  parameter  R  is  the  Reynolds  number  which  is  a  non-dimensional  measure  of 
the  ratio  of  inertial  forces  to  viscous  forces.  The  equatlona  (1.1)  hold  for  t  in 
(0,T)  and  x  in  a  domain  Q  in  r",  and  on  the  boundary  of  Q  there  are  n  acalar 
boundary  conditions.  For  definiteness  we  will  consider  the  boundary  conditions  to  be 
Dirlchlet  conditions  on  the  velocity,  i.e. 

(1.2)  u(x,t)  “  b(x,t)  on  3Qx[0,T] 

Much  of  what  we  say  will  not  depend  on  the  form  of  the  boundary  condition. 

There  are  a  number  of  systems  similar  to  the  Navier-stokes  equations  (1.1)  idilch  are 
aleo  of  significant  Interest.  Among  these  are  the  Bossineeq  equations  describing  fluids 
with^small  temperature  gradients  and  buoyancy  effects,  equations  of  non-Newtonian  fluid 
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flow,  and  tha  Stokes  aquations  describing  slow  viscous  flow,  Bte  Stokes  aquations  Bay  be 
written  as 


(1.3) 


±  „2* 

-  Vp  -  V  u 

MS  -  0  . 


in  a  x  [0,  T] . 


For  the  mathematical  theory  of  the  systeaa  (1.1)  and  (1.3)  we  refer  to  Ladyshenskaya 
(1963)  and  Teoaa  (1979). 


We  begin  our  discussion  of  numerical  methods  by  considering  what  features  are 
required  in  a  solution  algorithm.  First,  nusMrical  solutions  are  often  required  for  non¬ 
simple  gemetries  which  involve  curvilinear  grids  arising  from  boundary-fitted  coordinate 
systems.  Also,  the  grids  My  be  stretched  to  resolve  boundary  layers  or  other  flow 
phenomna  of  interest.  Secondly,  both  the  velocity  and  pressure  should  be  accurately 
determined.  The  pressure  can  be  mathematically  eliminated  from  (1.1)  and  (1.3)  by 
projecting  the  first  equation  in  each  system  onto  the  apace  of  divergence  free  vector 
fields.  However  the  determination  of  the  pressure  is  very  important  for  physical 
applications,  in  fact,  it  is  often  the  primary  objective.  One  is  frequently  interested  in 
the  pressure  force  induced  on  an  object  by  a  given  velocity  distribution. 

These  considerations  show  that  what  is  desired  for  the  Incompressible  Navler-Stokes 
equations  are  schemes  which  give  accurate  determination  of  the  velocity  and  pressure  for 
non- recti linear  grids.  It  is  also  desirable  to  have  methods  that  easily  extend  to  more 
general  systems  such  as  the  Boasinesq  equations. 

A  great  many  schemes  have  been  proposed  for  solving  the  incompressible  Navier-Stokes 
equations  and,  in  a  review  such  as  this,  it  is  impossible  to  acknowledge  more  than  a  small 
fraction  of  them.  We  have  tried  to  consider  only  questions  of  a  more  generel  nature 
without  dwelling  on  particular  features  of  the  scheMS.  Thus  we  refer  to  only  those 
papers  which  were  the  first  to  introduce  new  Mthods  or  ideas.  In  particular,  the  many 
fascinating  papers  dealing  with  applications  of  scheaas  to  probleM  of  scientific  interest 
are  not  mentioned  at  all. 


1 
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N*  apologise  to  thoao  4iom  work  we  htv*  overlooked.  Any  notification  of  ouch 
overalght  will  be  appreciated  and  will  be  corrected  In  a  eore  coeplete  review  which  the 
author  lntenda  to  write  later.  Me  hope  nonatheleae  that  this  review  will  etiaulate 
further  work  on  this  important  topic. 
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2.  Modifications  of  tha  Naviar-Stokaa  Equations 

The  Naviar-Stokaa  aquations  (1.1)  ara  frequently  modified  with  tha  intent  of 
obtaining  a  for*  which  is  aore  amenable  to  numerical  solution.  In  this  section  we  examine 
several  of  these  modifications. 

One  common  way  of  modifying  the  Navier-Stokea  equations 

(1.1)  is  to  replace  it  by  the  system 

♦  (i  •  $)u  ♦  “  V2u  -  f(x,t) 

(2.1) 

9*p  -  M  -  l  u1  u3 
i.j  Xj  Xi 

The  last  equation  of  (2.1)  ia  obtained  by  taking  the  divergence  of  the  first  equation  of 

(1.1)  and  using  the  last  equation  of  (1.1)  to  eliminate  the  divergence  of  the  velocity. 

The  motivation  for  using  the  system  (2.1)  is  the  hope  that,  when  discretised,  the  pressure 
can  be  obtained  using  standard  methods  for  inverting  the  diecrete  Laplacian.  However  the 
eystem  (2.1)  haa  a  grave  disadvantage  in  that  it  requires  an  additional  boundary 
condition.  Several  suggestions  have  been  made  for  this  boundary  condition. 

One  suggestion  is  that  the  additional  boundary  condition  be  given  by  the  normal 
derivative  of  the  pressure  as  determined  by  the  first  equation  of  (1.1)  or  (2.1)  evaluated 
on  the  boundary,  e.g.  Roache  (1972,p  194).  The  formula  is 

(2.2)  "  n**ut  +  ^  v*u  -  f(x,t)). 

This  however  is  not  satisfactory  as  a  boundary  condition  since  it  is  not  independent  of 
the  system  of  differential  equations  (2.1),  and  thus  leaves  the  overall  system 
underdetermined . 

Another  boundary  condition  triiich  is  used  along  boundaries  corresponding  to  physical 
surfaces  is  to  set  the  normal  derivative  of  the  pressure  equal  to  zero,  which  is  valid  in 
the  limit  for  high  Reynolds  number  flows.  This  gives  a  well-determined  system,  however, 
its  solutions  are  not  in  general  solutions  of  (1.1). 


X  good  My  to  illuatrata  thaaa  ideaa  la  to  oonaldar  ateady  tvo-dlmenaional  flow  In  a 


ehannal.  In  tha  region  0<x<L,  0<y<1  we  have  the  exact  aolution  to  (1.1) 

u  ■  y(1-y) 

(2.3)  . 

v  ■  0 


whore  u  -  (u,v)  la  the  velocity.  We  regard  thla  eolation  aa  determined  by  the 
apecifleatlon  of  the  velocity  oomponenta  on  tha  boundary,  l.e.  (1.2).  Dalng  the  preaaure- 
equation  ayatem  (2.1),  If  one  apacifiae  the  normal  derivative  of  preeaure  then  one  neede 
to  know  the  exact  value  of  dp/dx  or  one  obtaina  a  aolution  to  (2.1)  which  la  different 
from  (2.3)  and  la  not  a  aolution  of  (1.1). 

On  the  other  hand,  the  boundary  condition  (2.2)  underapedfiea  the  aolution  aince  any 
aolution  of  (2.1)  with 

S’  ♦<*'*> 

for  arbitrary  4  aatiafiea  (2.1)  and  (2.2). 

Thla  elmple  example  ehowo  that  aolving  (1.1)  via  (2.1)  la  unlikely  to  give  correct 
aolutiona.  Indeed  the  mathoda  uaing  (2.1)  or  aiallar  ayatama  have  difficulty  with 
accurate  preaaure  fielda  and  with  aatiafying  the  lncoaqpreealbility  condition  on  the 
velocltiaa  (aee  for  example  the  work  by  Boney,  Hefner,  Hirah,  and  Zoby  reported  in  Rubin 
and  Harrie  (1975)).  Roacha  (1972)  haa  a  diacuaaion  of  the  difficultlea  of  Obtaining  a 
aero  divergence  for  the  velocity  field  when  uaing  the  above  approach  (eee  also  Harlow  and 
welch  (1964)). 

The  only  boundary  condition  of  which  thla  author  knowa  which  will  make  (2.1) 
equivalent  to  (1.1)  ia  to  apecify  the  divergence  of  the  velocity  on  the  boundary.  Thla 
will  make  (2.1)  equivalent  to  (1.1)  for  low  Reynolda  numbera.  However,  thla  deatroye  the 
poaelbillty  of  iteratively  aolving  (2.1)  by  aolving  the  Polaaon  equation  for  preaaure 
aince  there  ia  no  boundary  condition  on  preaaure.  Becauae  of  theae  difficultlea  it  eeema 
beat  to  avoid  the  uae  of  (2.1)  and  eimilar  "pres aura- equation"  approachea  to  aolve  the 
Navlar-8tokea  equatione  (1.1). 
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A  number  of  researchers  make  transformations  similar  to  those  that  give  rise  to  (2.1) 


but  use  the  discrete  system  of  equations  obtained  from  either  a  finite  difference  or 
finite  element  approximation  of  (1.1).  It  is  usually  difficult  to  discern  enough  details 
from  their  writings  to  know  if  their  modified  equations  are  equivalent  to  their  original 
discrete  equations.  It  would  seem  that  most  of  these  discrete  pressure-equation  methods 
would  have  the  same  difficulties  as  the  differential  pressure-equation  approach,  that  is, 
the  number  of  possible  solutions  is  increased  and  an  additional  boundary  condition  is 
Introduced  to  guarantee  uniqueness.  Difficulties  would  arise  from  having  either 
inaccurate  boundary  conditions,  e.g.  normal  derivative  of  pressure  set  to  aero,  or  a 
boundary  condition  consistent  with  (2.2)  which  could  cause  ill-conditioning.  Chorin 
(1967)  discusses  other  shortcomings  of  this  approach.  In  any  case,  this  matter  is  more 
important  than  the  cursory  treatment  it  usually  receives  in  the  algorithm  descriptions. 

Another  approach  to  solving  the  steady  incompressible  Navier-Stokes  equations  is  the 
artificial  compressibility  method.  The  basic  idea  is  to  solve  a  time-dependent  system  of 
equations,  whose  steady-state  solutions  are  the  same  as  the  steady-state  solutions  of 
(1.1),  until  a  steady  state  is  reached.  Methods  have  been  proposed  by  Chorin  (1967), 
Yanenko  (1967),  and  Aubert  and  Deville  (1983).  The  convergence  rate  of  these  methods  is  a 
dependent  on  the  choice  of  time-discretization  used  to  solve  the  system.  A  difficulty 
which  arises  is  that  the  finite  difference  equations  may  not  have  a  steady-state 
solution.  The  reason  for  this  is  discussed  in  section  4. 

Another  common  method  is  to  use  the  'parabolized*  Navier-Stokes  equations  in  which 
the  second-derivatives  in  the  stream-wise  direction  are  removed.  Because  of  its  limited 
applicability  and  uncertain  justification  we  will  not  discuss  this  method  here  except  to 
note  that  often  an  analogue  of  (2.1)  is  derived  and  thus  some  of  our  observations  on  (2.1) 
also  apply  to  the  parabolized  equations.  Baithby  and  Schneider  (1979)  discuss  these 
difficulties  for  three-dimensional  flow  problems. 

A  popular  method  for  two-dimensional  and  axisymmetric  flows  is  the  stream- function 
and  vorticity  formulation.  For  two  dimensions,  the  stream  function  4>(x,y,t)  is  defined 
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IWo  boundary  condition*  ara  required  by  this  ayatan  and  corresponding  to  (1.2)  thay  ara 

given  by  specifying  both  components  of  tha  gradiant  of  4*  Note  however  that  tha  ayataai 

(2.4)  would  be  anich  easier  to  solve  if  one  had  one  boundary  condition  on  4  and  one  on 

(.  In  fact,  what  ia  usually  dona  in  finite  difference  calculations  is  that  (2.4b)  is 

2 

used  to  supply  values  of  (  on  tha  boundary,  and  tha  second  noraal  derivative  in  7  f  is 
approximated  using  the  boundary  conditions  for  tha  gradient.  The  nonlinear  system  is  then 
solved  iteratively,  (2.4a)  deteraing  C  given  4  and  (2.4b)  determining  4  given  (■ 
This  approach  works,  but  the  boundary  condition  on  C  causes  great  difficulty.  That  this 
oust  be  so  is  easy  to  see.  Suppose  at  soae  iteration  there  are  errors  in  4  of  aagnitude 

e.  Then  using  a  discrete  version  of  (2.4b)  as  the  boundary  condition  for  C ,  the  next 

-2 

iterate  of  c  will  have  errors  of  magnitude  c  h  where  h  ia  the  noraal  grid 
_2 

spacing.  (The  h  cones  fraa  the  second  noraal  derivative  in  the  laplacian.)  Thus  the 

errors  in  4  become  magnified  in  C .  The  usual  way  to  treat  this  is  to  under- relax  the 
boundary  values  of  C  on  the  boundary.  That  is 
cn+’  ■  (1-U)Cn  ♦  U  (-vy*1)  , 

2 

where  as  we  have  seen  u  nust  be  proportional  to  h  ,  giving  a  very  slow  rat*  of 
convergence.  In  spite  of  being  ao  slow  this  procedur*  la  often  need  and  can  give  accurate 
solutions. 

If  on*  eliminates  the  vorticity  from  the  system  (2.4)  on*  obtains  a  single  aquation 
which  for  steady  flow  is  a  fourth  ordar  elliptic  equation  in  4, 
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(2.5) 


0 


;  tS  -  ♦  v2*  +  *  v2*  - 

R  y  x  x  Ty 

This  equation  requires  two  boundary  conditions  and,  equivalent  to  (1.2),  they  are  given  by 

specifying  t  and  its  normal  derivative.  This  gives  a  very  well-conditioned  numerical 
problem.  Recent  work  using  (2.5)  has  been  done  by  Schreiber  and  Xeller  (1983a).  Once  a 
solution  is  obtained  for  (2.5),  the  pressure  can  be  obtained  by  integrating  the  equation 
for  the  gradient  of  the  pressure  from  (1.1).  If  the  equation  (2.5)  is  used  with  non- 
uniform  grids  and  coordinate  transformations,  the  discretised  system  can  be  quite 
formidable. 

Based  on  the  above  observations  it  would  seem  that  of  the  methods  we've  discussed 
only  the  original  version  of  the  Navier-Stokes  equations  (i.e.  (1.1))  and  the  stream 
function  formulation  (2.5)  can  lead  to  well-conditioned  numerical  methods. 

3.  Finite  Difference  Schemes 

In  this  section  we  discuss  finite  difference  schemes  for  the  Navier-Stokes  equations 
(1.1)  in  the  form  given,  we  will  concentrate  on  schemes  which  are  formally  second-order 
accurate,  although  there  are  a  number  of  schemes  triilch  are  first-order  accurate  (e.g. 
Kzivickil  and  Ladyzhenskaya  (1966),  Temam  (1979)). 

We  consider  first  the  so-called  staggered  mesh  schemes.  The  second-order  accurate 
staggered  mesh  scheme  for  a  uniform  cartesian  grid  assigns  the  values  of  each  of  the 
velocity  components  and  the  pressure  to  different  interlaced  grids.  In  two  dimensions 
with  velocity  components  u  and  v,  one  may  assign  values  of  u  to  grid  locations 
((1  +  */£  )h,jh),  values  of  v  to  (ih,(j  )h),  and  values  of  p  to  (ih,jh),  e.g. 
Harlow  and  Welch  (1965),  Patankar  and  Spalding  (1972),  Raithby  and  Schneider  (1979), 

Brandt  and  Dinar  (1979).  This  method  works  very  well  as  long  as  the  geometry  is 
rectangular  and  the  grid  is  uniform.  Non-uniform  grids  and  grid  supping  techniques  cannot 
be  conveniently  handled,  although  Liu  and  Krause  (1979)  have  developed  a  staggered  mesh 
scheme  for  use  with  general  geometries.  The  staggered  mesh  schesws  also  have  some 
difficulty  at  boundaries.  For  example,  when  both  velocity  components  are  specified  at  a 
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boundary  than  that  velocity  component  whose  mesh  lines  do  not  lie  on  the  boundary  requires 
some  special  treatment .  A  method  for  treating  curved  boundaries,  has  been  developed  by 
Viecelli  (1971). 

The  central  difference  scheme  on  a  uniform  rectangular  mesh  assigns  values  of  all  the 
variables  to  each  grid  point*  The  divergence  and  gradient  operators  are  approximated 
using  central  differences  and  the  laplacian  is  approximated  by  the  standard  five-point 
discrete  Implacian.  Central  difference  schemes  have  been  introduced  by  Chorin  (1967, 

1968)  for  time-dependent  calculations,  and  Aubert  and  Deville  (1983)  have  developed  a 
fourth  order  accurate  central  difference  schesm  for  steady  flows. 

An  important  concept  for  finite  difference  schemes  for  elliptic  systems  such  as  (1.1) 
and  (1.2)  is  that  of  regularity  (see  Bube  and  Strikwerda  (1983),  and  also  Prank  (1968), 
Brandt  and  Dinar  (1979)).  Regular  schemes  give  rise  to  regularity  estimates  analogous  to 
those  in  the  theory  of  elliptic  systems  of  differential  equations.  Solutions  to  regular 
difference  schemes  will  in  general  be  smoother  than  solutions  to  non-regular  schemes  and 
also  will  be  more  accurate  approximations  to  the  solutions  of  the  differential 
equations.  The  central  difference  scheme  is  non-regular  (Bube  and  Strikwerda  (1983)),. 
which  results  in  non-smooth  solutions.  The  lack  of  smoothness  is  most  noticeable  in  the 
pressure.  The  staggered  mesh  scheme  however  is  regular.  The  advantage  of  the  central 
difference  scheme  is  that  it  is  easily  implemented  with  non-uniform  grids  as  introduced  by 
coordinate  changes. 

It  should  be  emphasised  that  none  of  the  difficulties  mentioned  above  are 
insurmountable.  Both  the  staggered  mesh  and  central  differencing  schemes  have  been  used 
and  often  quite  successfully. 

A  scheme  which  incorporates  both  the  ease  of  use  which  central  differences  possess 
and  regularity  is  the  regularized  central  difference  scheme  which  has  been  introduced  by 
the  author  (Strikwerda  (1983)).  This  scheme  has  been  shown  to  be  second-order  accurate 
for  both  the  velocity  and  pressure  when  applied  to  the  Stokes  equations  (1.3)  (Strikwerda 
(1983)).  In  this  scheme  the  derivatives  of  pressure  are  approximated  as 
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3p  ,  2,  .2 

3^  ■  *kOP  '  0hk6k-\>P 


and  the  first  derivatives  of  the  velocity  in  the  divergence  equation  are  approximated  as 


3u  .  k  ,2.  ,2  k 

r —  “  6.  .u  “  oh.  o.  .5.  u  , 
3x^  kO  k  k+  k- 


where  a  is  a  non-zero  constant  and  .  and  6  are  the  centered,  forward,  and 

ku  k*  k- 

backward  divided  differences,  respectively.  The  laplacian  is  approximated  with  the  usual 
five-point  stencil.  Note  that  for  a  equal  to  one-sixth  the  approximations  (3.1)  and 

(3.2)  are  third-order  accurate. 

Since  the  regularized  central  difference  scheme  is  a  variant  of  the  central 
difference  scheme  it  is  easy  to  implement  with  coordinate  changes  and  boundary-fitted 


4.  The  Integrabilitv  Condition 

In  order  for  a  solution  to  exist  to  the  system  (1.1)  with  boundary  conditions  (1.2) 
it  is  necessary  that  the  boundary  data  satisfy  the  integrability  condition.  This 
condition  is  obtained  by  the  divergence  theorem  as 


(4.1)  0  =  /  /  $*u  ■  /  n*u  “  /  n*b  . 

“  3(1  3(1 


Thus  if  the  Integral  of  the  normal  component  of  the  boundary  data  is  non-zero,  there  is  no 
solution. 

Now  then  one  discretizes  the  equations  (1.1)  one  obtains  a  discrete  analog  of  the 
incompressibility  condition 

(4.2)  Dn*uh  “  0  in  (lh 
and  the  boundary  condition 

(4.3)  Uh  -  bn  on  3(lh. 


( 


Thera  is  aleo  a  discrete  integrability  condition  analogue  to  (4. 1 )  which  oust  be 
satisfied. 

There  are  at  least  two  ways  to  satisfy  the  discrete  intagrability  condition,  the 
first  method  would  be  to  analyse  the  aatrlx  corresponding  to  (4.2)  and  determine  the  nail 
space  of  the  adjoint  matrix.  If  the  data  is  constrained  to  be  orthogonal  to  this  nail 
space  then  a  solution  will  exist.  This  approach  is  impractical  for  many  situations , 
especially  if  coordinate  changes  have  been  employed  since  then  the  matricee  are  not  easy 
to  analyse.  However,  Amit,  Hall,  and  Forshlng{1982>  have  given  a  method  which  determines 
the  null  space  of  the  discrete  divergence  operator  for  special  rectilinear  gride  for  two- 
dimensional  problems. 

A  second  approach,  introduced  by  Strikwerda  (1983),  is  to  replace  (4.2)  by 

(4.4-,  V%  -  ah 

where  6.  is  a  constant  chosen  to  guarantee  a  solution.  The  value  of  6.  must  be 
n  n 

determined  as  part  of  the  solution.  As  shown  in  the  examples  of  Strikwerda  (1983)  i.  is 

n 

at  least  0(h2)  for  the  regularised  central  scheme. 

Another  way  of  looking  at  the  condition  (4.4)  is  as  follows.  As  shown  by  TOmem 

(1979)  and  others,  any  discrete  divergence  operator  D.  ,  defined  only  on  the  interior  of 

n 

the  grid,  has  a  corresponding  gradient  operator  6^  defined  by 

(4.5)  <D  •£,♦)  ♦  (u,Gi»)  -  0 

n  n 

♦ 

for  all  grid  vector  function  u  end  scalar  functions  9  which  vanish  on  the  boundary. 

If  one  wishes  to  satisfy  (4.2)  at  each  point  of  the  interior  then  (4.5)  with  9  taken  to 
be  one  at  each  Interior  point  gives  the  requirement  that 

(4.6)  (D  *u, 1 )  +  (u.Gil)  -  0 

n  n 

must  be  satisfied.  This  formula  is  the  analog  of  the  integrability  condition  (4.1),  and 

♦ 

the  second  term  in  (4.6)  will  usually  Involve  only  the  values  of  u  on  and  near  the 
boundary.  If  the  constraint  (4.6)  is  not  satisfied  then  the  data  must  be  modified  so  that 

(4.6)  is  satisfied.  The  use  of  (4.4)  in  place  of  (4.2)  is  one  say  by  irtiich  (4.6)  can  be 
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satisfied.  An  advantage  of  using  (4.4)  over  approaches  which  would  Modify  the  boundary 
data  of  u  is  that  (4.4)  requires  no  explicit  knowledge  of 

Xt  is  interesting  to  note  that  for  the  staggered  aesh  scheae  on  a  unifora  grid  one 
can  easily  satisfy  the  discrete  integrability  condition  since  the  calculus  of  finite 
differences  aiaics  the  differential  calculus  very  closely,  see  e.g.  Kzivickii  and 
Ladyzhenskaya  (1966).  Similar ly,  Chorin  (1969)  proves  the  convergence  of  a  central 
difference  aches*  for  the  time-dependent  Navier-Stokes  equations  on  a  periodic  rectangular 
sash.  An  essential  eleaent  of  these  proofs  is  that  one  has  a  convenient  fora  of  the 
finite  difference  analogue  of  the  divergence  theorea  of  the  differential  calculus.  Liu 
and  Krause  (1979)  develop  a  staggered  grid  scheae  for  non-rectangular  grids  and  the 
success  of  their  scheae  is  due  to  their  careful  treataent  of  the  integrability  constraint. 

S.  Solving  the  Finite  Difference  aquations. 

Having  deterained  which  fora  of  the  Navier-Stokes  equations  is  to  be  used  and  given 
the  finite  difference  equation*  to  be  eaployed,  one  is  still  left  with  finding  an 
efficient  aethod  for  solving  the  equations.  We  have  already  discussed  solving  the 
vortlcity-streaaf unction  formulation  in  section  2.  The  fonaulation  using  just  the 
streaafunction  can  be  solved  by  iterative  aethods  which  employ  direct  solvers  for  the 
biharaonic  equation  or  by  Newton-like  methods  (Schreiber  and  Keller  (1983a)). 

For  the  primitive  variable  formulation  the  basic  problea  reduces  to  solving  systems 
of  the  fora 

A(u)  u  ♦  <Sp  ■  b, 

(5.1)  Du  -  bj 

where  A(u)u  arises  froa  discretizing  the  velocity  portion  of  the  first  equation  of  (1.1) 
and  6  and  D  are  the  discrete  gradient  and  divergence  respectively.  A  coaaon  aethod  to 
solve  (5.1)  is  to  use  an  iterative  method  each  iteration  consisting  of  one  step  of  an 
Iterative  swthod  appropriate  for 
A(u)  w  •  ct 


wt 
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-1 


together  with  a  pniiuri  update  of  tha  for* 


wl 


v  _  v+1 
p  -  YDu 


Savaral  methods  hava  boon  proposed.  The  first  seems  to  have  boon  given  by  Chorin  (1968) 
for  tha  t is •- dependent  Mavis r-8tokes.  Other  Methods  have  been  introduced  for  finite 
elssMnt  Methods,  see  e.g.  Fortin  end  Glowinski  (1982),  end  by  Strikwerde  (1983)  for  finite 
difference  schemes. 

Recently  an  intsresting  Method  for  soXving  the  finite  difference  equations  arising 
from  the  Mavis r-8t okas  aquations  has  been  advanced  by  tait,  HaXX ,  and  Foraching 
(1982), (see  aXao  Ragaaan  (1982)).  An  easentiaX  requirement  is  that  one  be  able  to 
determine  the  null  space  of  tha  discrets  divergence  operator.  If  this  is  possible  then 
their  method  appears  to  work  wall.  At  presant,  this  Method  seems  limited  to  certain 
discretisations  of  the  rectangle. 


6.  Determining  the  Accuracy  of  a  Method. 

With  these  various  methods  available  an  obvious  question  arises  as  to  how  one  decides 
which  methods  are  more  accurate  or  suitable  than  others.  This  is  not  an  easy  question  to 
answer. 

Most  of  the  numerical  methods  tdiich  have  been  presented  have  been  tested  on  the 
'driven  cavity  problem.  However,  until  recently  there  has  not  been  a  numerical  solution  of 
this  problem  tdiich  is  of  sufficiently  high  accuracy  to  be  a  reliable  comparison.  The  work 
of  8chrelber  and  Keller  (1983a)  should  provide  such  a  standard  by  which  accuracy  can  be 
determined  for  other  methods.  Tho  experimental  data  on  this  problem  does  not  appear 
adequate  to  provide  good  checks  of  accuracy.  Xt  should  also  be  mentioned  that  seme 
methods  which  were  designed  to  do  well  on  flows  like  the  driven  cavity  give  totally 
Incorrect  solutions  due  to  excessive  artificial  viscosity  (Strikwerde  (1982)),  wee  also 


Schreiber  and  Keller  (1983b) 


A  teet  problem  with  an  analytic  solution  to  tha  time-dependent  Navier-Stokes 
aquations  waa  uaad  by  Chorln  (1967).  This  problem  could  ba  usad  aora  widaly  aa  a  taat 

casa. 

Tha  advantage  of  using  analytic  solutions  to  taat  numerical  schemes  la  that  this  is 
the  only  way  one  can  determine  tha  true  accuracy  of  tha  scheme.  As  shown  in  section  2 
schemes  which  use  tha  pressure  aquation  approach  can  ba  reroth-order  accurate,  similarly 
for  certain  upwind  differencing  schemes  (sea  Strlkwerda  (1982)).  If  these  scheaMS  had 
bean  adequately  tested  either  they  would  not  have  been  published  or  we  would  know  how  such 
confidence  to  put  in  them. 

It  is  very  easy  to  obtain  analytic  solutions  to  the  system  (1.1).  For  example,  for 

any  aauoth  function  9(x,y,t)  one  can  generate  solutions  in  two  dimensions  by  prescribing 

(6.1)  u  -  *  ,  v  -  -* 

y  * 

and  setting  the  pressure  arbitrarily.  The  force  f(x,y,t)  is  determined  so  as  to  make 
(6.1)  a  solution.  For  example,  with  (6.1)  one  can  sat 
P  <*x2  ♦  ♦/> 
f,  -  -  ;(V28)y  - 

f ,  -  ♦  ^(V2*)  -  ♦  9  V2* 

2  r  x  y  y 

or 

p  -  'h  <♦  2  ♦  *  2>  ♦ 

*  y 

*,  -  -  i«v2+>y  -  *«’2*)x 

f2  •  ♦  -  t<72*>y» 

A  paper  idiich  is  notable  for  its  careful  analysis  of  the  accuracy  is  that  by  Aubert 
and  Devllle  (1983).  They  present  a  acheam  which  is  fourth-order  accurate  in  the 
velocities  and  third-order  in  the  pressure.  The  accuracy  is  well  demonstrated  with  a 
series  of  test  problems. 
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7.  Conclusion* 

In  this  review  m  have  examined  method  a  of  solving  ths  incompressible  Nevier-S token 
aquations.  It  is  hoped  that  ths  consents  and  observations  given  hers  sill  stimulate 
others  to  develop  sore  accurate  and  efficient  nunarical  aathoda  for  these  important 
equations . 
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